---
title: "Regression"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output: html_document
---

```{r setup, include=FALSE}
rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
source("functions.R")
require(memisc)
require(stringi)
library(stargazer)
library(coefplot)
library(interplot)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(stargazer)
library(stringr)
library(tidyverse)
library(ggplot2)
library(psych)
library(margins)
library(dplyr)
library(modelsummary)
library(makedummies)
library(ggpubr)

theme_set(theme_sjplot())

# regression setting
dat <- cbind(readRDS("data/data_lss.RDS"), readRDS("data/data_dictionary.RDS"))

events <- yaml::read_yaml("events.yml")

# create hostile relationship dummy
dat$hostile <- NA
dat$hostile[dat$date <= as.Date("2016-02-09")] <- "pre"
dat$hostile[as.Date("2016-02-10") <= dat$date] <- "post"
dat$hostile <- factor(dat$hostile, levels = c("pre", "post"))

# create Afghan-rule dummy
dat$Afghan_rule <- NA
dat$Afghan_rule[dat$date <= as.Date("2021-08-15")] <- "pre"
dat$Afghan_rule[as.Date("2021-08-16") <= dat$date] <- "ruler"
dat$Afghan_rule <- factor(dat$Afghan_rule, levels = c("pre", "ruler"))

colnames(dat)[6]<-"LSS"
```

# corelation
```{r}
stargazer(as.data.frame(dat), type = "text")

dsrp_statsc <- dat |> 
  dplyr::select(LSS, publication, hostile, Afghan_rule)
descr(dsrp_statsc)
summary(dsrp_statsc)
table(dsrp_statsc$publication)

hist(dat$LSS)

ggplot() +
  geom_histogram(aes(x = dat$LSS), color = "white", bins = 30)+
  labs(x = "LSS", y = "Frequency")
```

# regressions
```{r}
# H1: difference in organizations
reg01 <- lm(LSS ~ publication, data = dat)
mtable(reg01,
       summary.stats=c("sigma","R-squared","F","p","N"))

# H2: hostile
reg02 <- lm(LSS ~ publication + hostile, data = dat)
mtable(reg02,
       summary.stats=c("sigma","R-squared","F","p","N"))

reg03 <- lm(LSS ~ hostile * publication, data = dat)
mtable(reg03,
       summary.stats=c("sigma","R-squared","F","p","N"))

reg04 <- lm(LSS ~ hostile, 
            subset(dat, publication == "Sumud"))
mtable(reg04,
       summary.stats=c("sigma","R-squared","F","p","N"))

reg05 <- lm(LSS ~ hostile, 
            subset(dat, publication == "Nabaa"))
mtable(reg05,
       summary.stats=c("sigma","R-squared","F","p","N"))

# H3: Afghan rule
reg06 <- lm(LSS ~ publication + Afghan_rule, data = dat)
mtable(reg06,
       summary.stats=c("sigma","R-squared","F","p","N"))

reg07 <- lm(LSS ~ Afghan_rule * publication, data = dat)
mtable(reg07,
       summary.stats=c("sigma","R-squared","F","p","N"))

reg08 <- lm(LSS ~ Afghan_rule, 
            subset(dat, publication == "Sumud"))
mtable(reg08,
       summary.stats=c("sigma","R-squared","F","p","N"))

reg09 <- lm(LSS ~ hostile, 
            subset(dat, publication == "Nabaa"))
mtable(reg09,
       summary.stats=c("sigma","R-squared","F","p","N"))
```

# regression table
```{r}
reg_table<- mtable("M1: ornanization" = reg01, "M2: relationship" = reg03, "M3: Afghan rule" = reg07,
               summary.stats=c("adj. R-squared","p", "N"))
print(reg_table)
write_html(reg_table, file = "result.html")
```

# marginal effect
```{r}
# H1
p1 <- plot_model(reg01, type = "pred", terms = c("publication"), 
           axis.labels = "", title = "M1: Marginal effect of Islamist organizations")
print(p1)

# H2
p2 <- plot_model(reg03, type = "int", 
           axis.labels = "", title = "M2: Marginal effect of hostile relations")
print(p2)

# H4
p3 <- plot_model(reg07, type = "int", 
           axis.labels = "", title = "M3: Marginal effect of Afghan rule")
print(p3)

ggarrange(p1,
          NULL,
          p2, 
          p3 + theme(axis.title.y = element_blank()),
          nrow = 2, ncol = 2, align = "hv",
          heights = c(0.4, 0.6),
          common.legend = TRUE,
          legend = "right",
          hjust = 0)
```
